#include "HY1C_out.h"
#include <stdio.h>
#include <stdlib.h>
#include "l12_proto.h"

static int xcal_wave[NBANDS];
static int nband = 0;
static double pixwt = 1.0;
static int detfac = 1;

int xcal_band(int wave)
{
    int iband = -1;
    int iw;

    for (iw=0; iw<nband; iw++)
        if (xcal_wave[iw] == wave)
            iband=iw;

    return(iband);  
}

double *get_xcal(l1str *l1rec, int type, int wave)
{
    static double rvs [XTNBAND][XTNMSIDE][XTNDET][MAXPIX];
    static double m12 [XTNBAND][XTNMSIDE][XTNDET][MAXPIX];
    static double m13 [XTNBAND][XTNMSIDE][XTNDET][MAXPIX];

    int bandnum;
    
    if ((bandnum = xcal_band(wave)) < 0) {

        double rvstab[XTNORDER][XTNDET][XTNMSIDE][XTNTIME];
        double m12tab[XTNORDER][XTNDET][XTNMSIDE][XTNTIME];
        double m13tab[XTNORDER][XTNDET][XTNMSIDE][XTNTIME];
        int16  yrtab [XTNTIME];
        int16  dytab [XTNTIME];
        double sectab[XTNTIME];

        int16 year, day;
        int32_t sensorID = l1rec->sensorID;
        double rvsc1[4];
        double rvsc2[4];
        double m12c1[4];
        double m12c2[4];
        double m13c1[4];
        double m13c2[4];
        double ap, ap2, ap3, x1, x2;
        double sec, wt;
        int32_t ic, ip, ib, im, id, it, it1, it2;
        int32_t norder, ntime, ndet, nmside;

        char  file   [FILENAME_MAX] = "";
        char  name   [H4_MAX_NC_NAME]  = "";
        char  sdsname[H4_MAX_NC_NAME]  = "";
        int32 sd_id;
        int32 rank;
        int32 nt; 
        int32 dims[H4_MAX_VAR_DIMS]; 
        int32 nattrs;
        int32 sds_id; 
        int32 start[5] = {0,0,0,0,0}; 
        int32 end  [5] = {1,1,1,1,1};
        int32 status;

        bandnum = nband++;
        xcal_wave[bandnum] = wave;

        sprintf(file,"%s%s%d%s",l1rec->input->xcalfile,"_",wave,".hdf");

        HY1C_out("Loading XCAL rvs and polarization sensitivities from %s\n",file);

        sd_id = SDstart(file, DFACC_RDONLY);
        if (sd_id == FAIL){
            fHY1C_out(stderr,"-E- %s line %d: SDstart(%s, %d) failed.\n",
            __FILE__,__LINE__,file,DFACC_RDONLY);
            exit(1);
        }

        status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Coefficients"),&norder);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
                __FILE__,__LINE__,"Number of Coefficients",file);
            exit(1);
	}
        if (norder > XTNORDER) {
	  HY1C_out("-E- %s Line %d:  mismatch in norder\n",__FILE__,__LINE__);
            exit(1);
	}
	  
        status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Detectors"),&ndet);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
                __FILE__,__LINE__,"Number of Detectors",file);
            exit(1);
	}
        if (ndet > XTNDET) {
  	    HY1C_out("-E- %s Line %d:  mismatch in ndet\n",__FILE__,__LINE__);
            exit(1);
	}

        status = SDreadattr(sd_id,SDfindattr(sd_id,"Number of Mirror Sides"),&nmside);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
                __FILE__,__LINE__,"Number of Mirror Sides",file);
            exit(1);
	}
        if (nmside > XTNMSIDE) {
  	    HY1C_out("-E- %s Line %d:  mismatch in nmside\n",__FILE__,__LINE__);
            exit(1);
	}
	  
        status = SDreadattr(sd_id,SDfindattr(sd_id,"Length of Time Series ndates"),&ntime);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading global attribute %s from %s.\n",
                __FILE__,__LINE__,"Length of Time Series ndates",file);
            exit(1);
	}
        if (ntime > XTNTIME) {
  	    HY1C_out("-E- %s Line %d:  mismatch in ntime\n",__FILE__,__LINE__);
            exit(1);
	}
	  
        strcpy(sdsname,"year");
        start[0] = 0;
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) yrtab);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
                __FILE__,__LINE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"day");
        start[0] = 0;
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        status = SDgetinfo(sds_id, name, &rank, dims, &nt, &nattrs);
        status = SDreaddata(sds_id, start, NULL, dims, (VOIDP) dytab);
        if (status != 0) {
            HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
                __FILE__,__LINE__,sdsname,file);
            exit(1);
        } else {
            status = SDendaccess(sds_id);
        }

        strcpy(sdsname,"M11");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        for (ic=0; ic<norder; ic++) for (id=0; id<ndet; id++) for (im=0; im<nmside; im++) for (it=0; it<ntime; it++) {
            start[0]=ic; start[1]=id; start[2]=im; start[3]=it;
            status = SDreaddata(sds_id, start, NULL, end, (VOIDP) &rvstab[ic][id][im][it]);
            if (status != 0) {
                HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
                    __FILE__,__LINE__,sdsname,file);
                exit(1);
            } 
        }
        status = SDendaccess(sds_id);

        strcpy(sdsname,"m12");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        for (ic=0; ic<norder; ic++) for (id=0; id<ndet; id++) for (im=0; im<nmside; im++) for (it=0; it<ntime; it++) {
            start[0]=ic; start[1]=id; start[2]=im; start[3]=it;
            status = SDreaddata(sds_id, start, NULL, end, (VOIDP) &m12tab[ic][id][im][it]);
            if (status != 0) {
                HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
                    __FILE__,__LINE__,sdsname,file);
                exit(1);
            } 
        }
        status = SDendaccess(sds_id);

        strcpy(sdsname,"m13");
        sds_id = SDselect(sd_id, SDnametoindex(sd_id,sdsname));
        for (ic=0; ic<norder; ic++) for (id=0; id<ndet; id++) for (im=0; im<nmside; im++) for (it=0; it<ntime; it++) {
            start[0]=ic; start[1]=id; start[2]=im; start[3]=it;
            status = SDreaddata(sds_id, start, NULL, end, (VOIDP) &m13tab[ic][id][im][it]);
            if (status != 0) {
                HY1C_out("-E- %s Line %d:  Error reading SDS %s from %s.\n",
                    __FILE__,__LINE__,sdsname,file);
                exit(1);
            } 
        }
        status = SDendaccess(sds_id);

        status = SDend(sd_id);


        for (it=0; it<ntime; it++) {
	    if (yrtab[it] > 0)
                sectab[it] = yds2unix(yrtab[it],dytab[it],0.0);
	    else 
	        break;
	}
        if (it != ntime) {
	  HY1C_out("-E- %s line %d: time field mismatch %d %d\n",__FILE__,__LINE__,it,ntime);
	    exit(1);
	}

        sec = yds2unix(*(l1rec->year),*(l1rec->day),*(l1rec->msec)/1000);
        for (it=0; it<ntime; it++) {
	  if (sectab[it] > sec) break;
	}
        if (it == 0)
	    it2 = it1 = it;
        else if (it == ntime) 
	    it2 = it1 = ntime-1;
        else {
            it2 = MAX(MIN(it,ntime-1),1);
            it1 = it2-1;
	}

        for (im=0; im<nmside; im++) for (id=0; id<ndet; id++) {

	    rvsc1[0] = rvstab[0][id][im][it1];
	    rvsc1[1] = rvstab[1][id][im][it1];
	    rvsc1[2] = rvstab[2][id][im][it1];
	    rvsc1[3] = rvstab[3][id][im][it1];
	    m12c1[0] = m12tab[0][id][im][it1];
	    m12c1[1] = m12tab[1][id][im][it1];
	    m12c1[2] = m12tab[2][id][im][it1];
	    m12c1[3] = m12tab[3][id][im][it1];
	    m13c1[0] = m13tab[0][id][im][it1];
	    m13c1[1] = m13tab[1][id][im][it1];
	    m13c1[2] = m13tab[2][id][im][it1];
	    m13c1[3] = m13tab[3][id][im][it1];

	    rvsc2[0] = rvstab[0][id][im][it2];
	    rvsc2[1] = rvstab[1][id][im][it2];
	    rvsc2[2] = rvstab[2][id][im][it2];
	    rvsc2[3] = rvstab[3][id][im][it2];
	    m12c2[0] = m12tab[0][id][im][it2];
	    m12c2[1] = m12tab[1][id][im][it2];
	    m12c2[2] = m12tab[2][id][im][it2];
	    m12c2[3] = m12tab[3][id][im][it2];
	    m13c2[0] = m13tab[0][id][im][it2];
	    m13c2[1] = m13tab[1][id][im][it2];
	    m13c2[2] = m13tab[2][id][im][it2];
	    m13c2[3] = m13tab[3][id][im][it2];

            if (it2 == it1) 
	        wt = 0.0; 
            else
                wt = (sec-sectab[it1])/(sectab[it2]-sectab[it1]);

            if (sensorID == HMODISA || sensorID == HMODIST) {
	      switch (l1rec->input->resolution) {
                case 250:
  	          pixwt  = 0.25;
                  detfac = 4;
                  break;          
  	        case 500:
  	          pixwt  = 0.5;
                  detfac = 2;
                  break;          
	      }
  	    } 

            for (ip=0; ip<l1rec->npix; ip++) {

	        ap  = (double) l1rec->pixnum[ip]*pixwt;
                ap2 = ap*ap;
                ap3 = ap2*ap;

                x1 = rvsc1[0]+ap*rvsc1[1]+ap2*rvsc1[2]+ap3*rvsc1[3];
                x2 = rvsc2[0]+ap*rvsc2[1]+ap2*rvsc2[2]+ap3*rvsc2[3];
                rvs[bandnum][im][id][ip] = x1 + wt*(x2-x1);

                x1 = m12c1[0]+ap*m12c1[1]+ap2*m12c1[2]+ap3*m12c1[3];
                x2 = m12c2[0]+ap*m12c2[1]+ap2*m12c2[2]+ap3*m12c2[3];
                m12[bandnum][im][id][ip] = x1 + wt*(x2-x1);

                x1 = m13c1[0]+ap*m13c1[1]+ap2*m13c1[2]+ap3*m13c1[3];
                x2 = m13c2[0]+ap*m13c2[1]+ap2*m13c2[2]+ap3*m13c2[3];
                m13[bandnum][im][id][ip] = x1 + wt*(x2-x1);
	    }   
	}
        
    }

    switch (type) {
      case XRVS: return (&rvs[bandnum][l1rec->mside][l1rec->detnum/detfac][0]); break;
      case XM12: return (&m12[bandnum][l1rec->mside][l1rec->detnum/detfac][0]); break;
      case XM13: return (&m13[bandnum][l1rec->mside][l1rec->detnum/detfac][0]); break;
      default: 
        HY1C_out("-E- %s line %d: bad xcal type specified\n",__FILE__,__LINE__);
        exit(1);
    }   

}
